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Dynamic particle-scale numerical simulations are used to show that the shear thickening observed 
in dense colloidal, or Brownian, suspensions is of a similar nature to that observed in non-colloidal 
suspensions, i.e., a stress-induced transition from a flow of lubricated near-contacting particles to a 
flow of a frictionally contacting network of particles. Abrupt (or discontinuous) shear thickening is 
found to be a geometric rather than hydrodynamic phenomenon; it stems from the strong sensitivity 
of the jamming volume fraction to the nature of contact forces between suspended particles. The 
thickening obtained in a colloidal suspension of purely hard frictional spheres is qualitatively similar 
to experimental observations. However, the agreement cannot be made quantitative with only 
hydrodynamics, frictional contacts and Brownian forces. Therefore the role of a short-range repulsive 
potential mimicking the stabilization of actual suspensions on the thickening is studied. The effects 
of Brownian and repulsive forces on the onset stress can be combined in an additive manner. The 
simulations including Brownian and stabilizing forces show excellent agreement with experimental 
data for the viscosity rj and the second normal stress difference A 2 . 


Introduction The rheology of dense suspensions is of 
considerable theoretical and technological importance, 
yet the shear rheology of even the simplest case of a sus¬ 
pension of hard spheres in a Newtonian suspending fluid 
is incompletely understood [T] . Many of the features ob¬ 
served in these suspensions, including shear thinning [5] 
or thickening mil] and the magnitudes and even the al¬ 
gebraic signs of normal stress differences [5], are at best 
understood at a qualitative level, and a general theoret¬ 
ical framework is lacking. Furthermore, there has been 
a tendency to treat the rheology of Brownian (colloidal) 
suspensions and non-Brownian suspensions as distinct. 

Recently, a picture has emerged in which central as¬ 
pects of the rheology of non-Brownian dense suspensions 
are interpreted as manifestations of proximity to jam¬ 
ming transitions in the parameter space. These transi¬ 
tions are singularities whose locations in the volume frac¬ 
tion (p and shear stress a plane depend on the details of 
the microscopic interactions (shape of the particles, fric¬ 
tion, interparticle forces). In turn, the locations of these 
singularities shape the large (p portion of the rheologi¬ 
cal landscape, i.e., the effective viscosity and the normal 
stresses as functions of (p and a. In particular, in the 
“stress-induced friction” scenario HEHn], shear thick¬ 
ening is a transition from a rheological response domi¬ 
nated by frictionless jamming to one controlled by fric¬ 
tional jamming upon increase of the shear stress. This 
transition is argued to be due to the creation of frictional 
contacts between particles at high stresses; the contacts 
are prevented at low stresses by a short-range stabiliz¬ 


ing repulsive force, as would be expected to be present 
to stabilize a colloidal dispersion against aggregation by, 
for example, attractive van der Waals forces [Ml dg. 
This picture contrasts with previous models for Brow¬ 
nian suspensions [TMT5] , in which frictional contacts are 
neglected based on idealized lubrication hydrodynamics. 

Suspensions are said to be colloidal, or Brownian, when 
the immersed particles are small enough: a commonly 
accepted upper bound for Brownian motion to be sig¬ 
nificant is a diameter of 1pm [2]. For these systems, 
the Brownian forces have been seen to be an essential 
factor in non-Newtonian behavior m- Physically, for 
a system of strictly hard colloidal spheres under shear 
in the Stokes regime, there are only two independent 
time scales: the inverse shear rate 7“^ and the diffu¬ 
sion time a^/Uo; here a is the sphere radius and Dg is 
the single-particle diffusion coefficient, which is related to 
the suspending fluid viscosity 770 and thermal energy A:bT 
through the Stokes-Einstein relation Dq = k^Tl&TTrjQa. 
The shear rate dependence of the rheology can be stated 
in terms of a competition between advection and diffu¬ 
sion described by the Peclet number Pe = dirrioa^j/kT. 
Smooth spheres with ideal lubrication resulting from hy¬ 
drodynamic interactions would be non-contacting, and 
hence exhibit the following rheology: a shear rate inde¬ 
pendent regime close to thermal equilibrium (that is for 
7 < T“^, where Tq is the “caging” time, or the typical 
time for which the thermal motion leads to a structural 
reorganization [20j ) where most forces are Brownian, fol¬ 
lowed by a shear thinning regime at intermediate values 
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of Pe, over which the Brownian forces become progres¬ 
sively less important relative to the hydrodynamic ones, 
and finally a shear thickening regime at large Pe that is 
dominated by the hydrodynamic lubrication forces due 
to increasingly smaller inter-particle gaps p^HTS] . di¬ 
verges at the glass transition (j) = (j)G, the system devel¬ 
ops a yield stress Uy above (/>g, and the low Pe viscosity 
plateau yields to an asymptotic shear thinning regime 
ry ~ CTy/y for 7 —>■ 0 pTM23j . Although numerical sim¬ 
ulations including only hydrodynamic interactions agree 
well with experimental data in the shear thinning regime, 
the simulated shear thickening regime is much weaker 
than is often experimentally observed, with the disagree¬ 
ment increasing as the volume fraction increases |18] . 

Most of the experimentally studied thickening suspen¬ 
sions are in the colloidal size range |1]. (Notable excep¬ 
tions include cornstarch suspensions.) It is thus essen¬ 
tial to address the validity of the stress-induced friction 
scenario for these systems. In this scenario, the num¬ 
ber of frictional contacts directly depends on the ratio of 
the shear stress to the Brownian stress scale aa^/ksT. 
At small stresses, i.e., when aa^/k^T <C 1, the thermal 
motion keeps particles separated and makes contacts un¬ 
likely. (In the equilibrium limit Pe —> 0, the average con¬ 
tact number per particle is zero at volume fractions be¬ 
low the jamming transition, otherwise the pressure would 
diverge as required by the virial equation for hard par¬ 
ticles [24].) When aa^/k^T 3> 1, however, the Brown¬ 
ian forces are not strong enough to overcome the forces 
bringing particles together due to the shear, and contacts 
are created. For dense suspensions, this shear activated 
friction mechanism can be related to a jamming tran¬ 
sition framework: at the largest shear stress for which 
shear thinning occurs, the rheology is dominated by the 
proximity to the frictionless jamming transition point at 
^j, in the sense that rheological properties are (roughly) 
diverging functions of the form (^ — ())j) “, whereas at 

shear stresses above shear thickening, the frictional jam¬ 
ming transition at dominates and rheological proper¬ 
ties scale with {p — with Aq and two positive 

exponents (whose exact values are still debated). Since 
0 j < 0j, this leads to shear thickening, which can be 
continuous or discontinuous depending on the proximity 
to <pj [m 125] . 

In this work, we show that simulations of frictional 
colloidal suspensions can reproduce both continuous 
and discontinuous shear thickening, hence demonstrat¬ 
ing that the stress-induced friction scenario extends to 
the Brownian case. Quantitative agreement with exper¬ 
iments cannot, however, be achieved with only frictional 
contacts and Brownian motion combined with hydrody¬ 
namic lubrication interactions. Instead, one must also 
consider the repulsive force that is induced between im¬ 
mersed colloids by the stabilization process. We study 
the qualitative influence of the range and amplitude of 



FIG. 1. The simulations consider a three-dimensional system 
of bidisperse spheres under simple shear with Lees-Edwards 
periodic boundary conditions. 

the repulsive force on the shear thickening. In particu¬ 
lar, we show that the effects of the Brownian and sta¬ 
bilizing forces on the onset stress are additive. For a 
suitable choice of amplitude and range of the repulsive 
force, the simulations of the relative viscosity and the 
second normal stress difference (which is large relative to 
the first) are in good agreement with recent experiments 
by Cwalina and Wagner pS] . 

Model description We study colloidal suspensions of 
hard spheres at positions X and velocities U in Stokes 
flow interacting through hydrodynamic (Fk), contact 
(Fc), stabilizing repulsive (^r), and Brownian (Fb) 
forces. We use a bidisperse system, with spheres of radii 
a and 1.4a at equal volume fractions. All our results are 
obtained with N = 500 particles in a cubic box with Lees- 
Edwards boundary conditions, as depicted in Fig.[T] The 
equation of motion is given by the following overdamped 
Langevin equation: 

0 = Fh -I- Fc -I- Fr -f Fb . (1) 

The hydrodynamic forces are the sum of a drag due 
to the motion relative to the surrounding fluid and 
a resistance to the deformation imposed by the flow: 
Fh = -Ffu(^) • {U - U°°) + RpEiX) : F°°, where 
17°° = 'jyiBx is the background imposed flow and E°° = 
{e^By + eye.x)'pl2 is the strain rate tensor. For the dense 
suspensions of interest in this work, it is assumed that 
the hydrodynamic interactions are dominated by near¬ 
contact lubrication interactions m and the long-range 
interactions are neglected. For smooth hard spheres, 
i?FU and Rfe contain lubrication terms that diverge 
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FIG. 2. Comparison of our simulation results (solid lines) for 
the viscosity as a fnnction of Pe with results obtained with 
Stokesian Dynamics by Foss and Brady [TS] (dashed lines) 
for a purely Brownian suspension at (j) = 0.45. We plot the 
total relative viscosity (black squares) as well as the individual 
hydrodynamic (blue triangles), Brownian (red diamonds) and 
contact (purple circles) contributions. 


with vanishing interparticle gap h. We keep only the 
leading order of these divergences (these are terms in h~^ 
and log/i“^, with h = h/a, for, respectively, the normal 
and tangential relative motion), and we regularize these 
terms with a microscopic cutoff <5 to mimic the parti¬ 
cle roughness [9]. [That is, the included terms scale as 
{h -I- 5/a) and log {h-\- 5/a) .] The results presented 
here are obtained with 5 = 10“^a. The results depend 
only weakly on the value of 5. The listing of the individ¬ 
ual matrix elements is given in El- 

Contacts are modeled by a linear spring consisting 
of both normal and tangential components, which is a 
simple model commonly used in granular physics El- 

Tangential and normal components of the contact force 

Hi) 

between two particles satisfy Coulomb’s friction 
law — I^|-^Cnor|- "^he Spring constants are cho¬ 

sen for each Pe and (j) so that the average minimal center- 
to-center distance dij between any two particles i and j 
is maintained as 1 — dij/{ai + dj) « 0.02 (here Ui and Uj 
are the particle radii). 

We take a stabilizing repulsive force that decays ex¬ 
ponentially with the interparticle gap h as I^’rI = 
F* exp(—/i/A), with a characteristic length A. This pro¬ 
vides a simple model of screened electrostatic interactions 
that can often be found in aqueous systems [55J [52] , in 
which case A is the Debye length. With one force scale 
and one length scale, this is the most basic parametriza- 
tion of a generic stabilizing force. 

The Brownian forces acting on different particles are 
correlated through the hydrodynamic interactions, and 
their second order cumulant is given by the fluctuation- 
dissipation theorem [30| . In a simulation with discretized 


time tn with intervals At, this becomes EIllSl] 

(FB(t™)FB(t„))-(FB(t„))(FB(t„)) = 

( 2 ) 

Owing to the dependence of the Brownian force term 
on the configuration as shown by Eq. the Langevin 
equation 0 contains a multiplicative noise. As a conse¬ 
quence, we must specify by which convention we evaluate 
the Brownian force |33j . A natural choice for numerical 
simulation is the Ito convention, which we adopt here; 
that is, when to = n in Eq. ([^, we evaluate i?FU at 
the beginning of the time step as RY\]{X(tn)). In this 
convention, the Brownian forces have a non-zero average 
(called “drift”) [321 E): 


(FB(t„)) = kBTRFv{X{tn)) - V - (X(t„)). (3) 

At each time step, we evaluate the contact and re¬ 
pulsive forces Fc and Fr, generate Brownian forces Fb 
according to Eqs. ([^ and ([^ (through the algorithm 
of |2I]), and we solve Eq. Q for the velocities: 

U-U^ =R^I-{Rbb-.E^+Fc + FbFFb). (4) 

Results We first show example results obtained for 
the “pure” Brownian case without repulsive forces, i.e., 
F* = 0. We show in Fig.[^ that the frictional contacts 
significantly modify the high Pe dynamics. In this fig¬ 
ure we compare the relative viscosity rji- at a volume 
fraction (j) = 0.45 obtained with our simulation using a 
friction coefficient /r = 1 with the result obtained using 
Stokesian Dynamics (with particles not making contact) 
by Foss and Brady [TH] as a function of Pe. There is 
very good agreement in the shear thinning regime at low 
Pe, where hydrodynamics and Brownian forces are dom¬ 
inant. At high Pe, however, where Stokesian Dynamics 
shows a very mild increase of the viscosity, our simula¬ 
tion shows a significant thickening that is due to contact 
stresses arising from the building up of a contact net¬ 
work. The hydrodynamic contribution in our simulation 
does not thicken at high Pe as predicted theoretically [36] 
for hydrodynamically-interacting hard spheres due to the 
lubrication cutoff, and instead becomes almost shear rate 
independent. 

The volume fraction dependence of the thickening due 
to the frictional contacts is shown in Fig.|^ We are able 
to obtain both continuous shear thickening for (j) ^ 0.55 
and discontinuous shear thickening for the highest vol¬ 
ume fractions </> > 0.55. The thickening is seen above 
an onset stress tTon ~ 5kBT/a^. This shows that the 
shear-induced friction mechanism introduced to explain 
the thickening of non-colloidal suspensions [8Hni[i3i ex¬ 
tends to the colloidal case, giving the same qualitative 
rheology. The values of the onset stress obtained exper¬ 
imentally are typically of order (Ton ~ lOQkBT/a^, how¬ 
ever I1S1I3III3HI, which is substantially larger than the 
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FIG. 3. Relative viscosity rj^ as a function of the Peclet number Pe (left) and as a function of the reduced shear stress (right) 
for a Brownian suspension of hard frictional spheres at several volume fractions 0.45 < 0 < 0.56. 




FIG. 4. Left: The effect of the repulsion amplitude F* on the rheology of a colloidal suspension for 0 = 0.5, for several 
values F* — 0, lO^feeT/a, 2 x 10^A:Br/a,5 x lO^ksT/a and W^ksT/a from light to dark color. The value of F* essentially 
influences the onset stress of thickening: increasing F* delays the thickening. As F* increases, the shear thinning regime is 
extended and the minimum value of the viscosity decreases slightly. Right: The effect of both the Brownian forces and the 
repulsive force can be summed up by plotting the relative viscosity as a function of the stress rescaled by the onset stress 
Uon = 5kBT/a^ + 0.01F’*/a^. This collapse of the data extends to the non-Brownian purely repulsive case (in black squares), 
for which tJon = O.OIF*/a^. 


value of about 5k^T/a^ found in our F* = 0 simulation. 
This can be understood as our simulation missing a repul¬ 
sive force that arises from the suspension “stabilization” 
in the experiments. 

We first look at the effect of the repulsive force on 
the thickening of the colloidal suspension. The relative 
viscosity curves for different values of F* are shown in 
the left panel of Fig.[^ for (p = 0.5. The main effect 
of the repulsion is, as expected, to push the onset of 
thickening to higher stresses. The relative viscosity in 
the thickened state is unaffected by the value of F*, as 
in this regime the repulsive force can be neglected relative 
to the hydrodynamic and contact forces. Note that the 
slope of the shear thinning is also the same for all of 
the simulated F*. In the right panel of Fig.|^ we show 
that the onset stress is approximately cJon ~ bk^T/a^ + 


Q.QIF*/a?. Thus, to a good approximation the effects of 
Brownian and repulsive forces on shear thickening can be 
combined in a simple additive manner. In this regard, the 
Brownian forces have an effect that is virtually identical 
to that of a potential repulsive force. 

Besides F*, our repulsive force contains another free 
parameter, the force decay length A. As we show in Fig.[^ 
A essentially controls (in conjunction with the Brownian 
motion) the strength of the shear thinning at low Pe. 
Increasing A, i.e., increasing the distance over which the 
repulsive force decays, makes the shear thinning more 
pronounced. This can be qualitatively understood: when 
A —> 0, the repulsive force disappears, and only the shear 
thinning due to the Brownian motion remains. 

We use simulations to assess the appropriate repulsive 
force to capture the behavior seen in experiment. We fo- 
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FIG. 5. The effect of the repulsion range A on the rheology 
of a colloidal suspension for (p = 0.5, from A/a = 0.01 to 
A/a = 0.05, from light to dark color. When A increases, the 
low shear rate viscosity increases, and as a consequence the 
shear thinning regime becomes more pronounced. The high 
shear rate rheology is unaffected by a change of A. 


cus here on the recent data by Cwalina and Wagner [26] , 
which include measurements of the shear stress and nor¬ 
mal stress differences for a suspension of silica beads with 
radius a = 260 nm in a low molecular weight polyethylene 
glycol (PEG) Newtonian suspending fluid ai T — 300 K. 
The particles are coated with octadecane chains to pro¬ 
vide steric stabilization [39]. The short range van der 
Waals attraction is also reduced by index matching be¬ 
tween particles and solvent. 

We obtained the best comparison with the results 
of Cwalina and Wagner [2S] by setting /r = 1, F* = 
5 X lO^keT/a and A = 0.02a, as shown in Fig.[^ The 
agreement with the experimental data for the relative 
viscosity is excellent. The second normal stress difference 
N 2 also shows a very good agreement with the experimen¬ 
tal data, being negative for all volume fractions; —N 2 I 0 
is in the range 0.15-0.4 for all Pe, as shown in Fig.[^ 
consistent with the behavior for non-Brownian suspen¬ 
sions m- Rather surprisingly, given the agreement for 
both 77r and N 2 , Ni disagrees substantively between our 
simulations and these experiments. In the experiments, 
Ni < 0, as predicted based on hydrodynamic force dom¬ 
inance |40j . Our simulations find weaker negative Ni 
for the lower volume fractions presented (</) = 0.50 and 
0.53), while for </> = 0.55, > 0. We note that 

Lootens et al. [41] observed a similar change in sign of Ni 
at the shear thickening transition. The parameters used 
in the simulations can be translated into SI units using 
the experimental parameters, with an inferred repulsive 
force at contact of F* « 79 pN and a repulsion range of 
A « 5.2 nm, which is to be compared to the thickness 
of the stabilizing polymer comb estimated to be 15 nm 
to 20 nm from the structure factor measured by neutron 
scattering |39j . The volume fractions used in the simu- 
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FIG. 6. Comparison with experimental data from Cwalina 
and Wagner [26] (black lines) for the relative shear viscosity 
(top), second (center) and first (bottom) normal stress differ¬ 
ence viscosities as functions of the Peclet number. Simulation 
results (colored lines and symbols) are obtained with a repul¬ 
sive force at contact F* = 5 x lO^kuT/a and a repulsion range 
A = 0.02a. 


lations to get the best agreement with experiments are 
always higher than the experimental ones. This might 
be attributed to the higher polydispersity used in the 
simulations, which lowers the viscosity at a given volume 
fraction. 
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FIG. 7. Second normal stress difference N 2 normalized by 
the shear stress as a function of the shear stress for several 
volume fractions 4> for the same conditions as in Fig.|^ 

Discussion Our simulation results suggest that the 
shear thickening of a dense suspension is fundamentally 
the same phenomenon for particles from 0(10nm) to 
O(lOOpm). Abrupt or discontinuous shear thickening 
occurs in suspensions close to jamming; as such, it is 
a geometric rather than a hydrodynamic phenomenon. 
(Hydrodynamics has been shown to provide a basis for 
weak continuous shear thickening in dilute |36j as well as 
concentrated nmn] suspensions.) It depends crucially 
on the existence of a mechanism preventing contacts at 
small stresses, such as Brownian motion in a purely col¬ 
loidal suspension. Interestingly, for understanding the 
qualitative rheological behavior, the Brownian force can 
be thought of as a potential repulsive force. A com¬ 
parison between our numerical results and experimental 
data [26j shows that in actual suspensions the repulsive 
effect of the Brownian force adds to an actual potential 
repulsive force. This analogy between Brownian force 
and repulsive force is not restricted to the thickening 
regime and has been proposed for suspensions at equilib¬ 
rium |42j and, recently, in the shear thinning regime [43]. 
Our work supports the view that a theoretical modeling 
of shear thickening should be centered on a framework 
common to Brownian and non-Brownian systems |44j . 
e.g., a geometric description [33111] including a stress- 
induced friction mechanism. 
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